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Fourier  sense.  This  paper  presents  a  general  purpose  numerical  treatment 
formulated  to  overcome  these  difficulties.  The  numerical  approach  is  based 
on  finite  difference  schemes  applied  in  conjunction  with  powerful  numerical 
ordinary  differential  equation  methods.  The  theory  is  examined  with  respect 
to  consistency,  stability,  and  convergence  of  these  numerical  procedures. 

A  numerical  example  is  included  to  demonstrate  the  validity  of  the  treatment 
Although  an  explicit  boundary  condition  is  absent  from  this  study,  a  derived 
boundary  condition  is  demonstrated  to  be  adequate. 
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A  NUMERICAL  TREATMENT  OF  THE  DYNAMIC  MOTION  OF  A 
ZERO  BENDING  RIGIDITY  CYLINDER  IN  A  VISCOUS  STREAM 

1.  INTRODUCTION 

Paidoussis*  worked  out  a  solution  to  the  dynamic  motion  problem. 

2 

Ortloff  and  Ives  studied  a  special  case  of  the  same  problem  and  expressed 
their  solution  in  the  form  of  an  infinite  series  involving  Gamma  and  Bessel 
functions.  Both  the  orders  and  the  arguments  of  Bessel  functions  are 
generally  complex  and  can  be  large  in  magnitude.  Furthermore,  evaluation  of  a 
Bessel  function  of  complex  order  is  laborious  and  time-consuming,  and  accuracy 
cannot  be  assured.  When  the  solution  proposed  by  Ortloff  and  Ives  is  applied 
to  the  nonhomogenous  problem  where  the  "upstream"  end  of  the  cylinder  is 
forced,  a  harmonic  time  dependence  is  assumed;  this  means  that  "forcing"  the 
system  by  an  arbitrary  time  function  will  require  multiple  solutions  combined 
in  the  Fourier  sense. 

To  overcome  these  difficulties,  a  general  purpose  numerical  approach  is 
introduced.  This  approach  discretizes  ,  n^r  ,  and  by  backward  and 

central  differences.  This  discretization  brings  the  dynamic  motion  equation 
into  a  system  of  second  order  ordinary  differential  equations.  This  system  is 
decomposed  into  a  system  of  first  order  ordinary  differential  equations.  A 
feasible  numerical  ordinary  differential  equation  method  is  then  used  to  solve 
this  system  with  optimal  efficiency. 

There  are  many  advantages  to  using  a  numerical  method  to  solve  the  problem 
of  dynamic  motion.  The  theory  is  well  developed  with  respect  to  consistency, 
stability,  and  convergence.  Numerical  methods  are  systematic  to  implement. 
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and  effective  techniques  can  be  used  to  accurately  accelerate  computations. 
When  a  numerical  approach  is  used,  the  laborious  evaluation  of  special 
functions  is  bypassed,  maximizing  accuracy  and  efficiency. 

This  report  begins  with  a  description  of  the  dynamic  motion  problem  and 
the  associated  initial  and  boundary  conditions.  A  numerical  approach  is 
introduced  and  the  supporting  theory  and  mathematical  formulation  are 
discussed.  An  example  is  given  to  demonstrate  the  validity  of  our  numerical 
solution  to  a  well  posed  dynamic  motion  problem.  The  computer  programs  are 
included  in  the  appendix. 
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2.  POSED  PROBLEMS 

The  motion  of  a  wire  suspended  in  a  fluid  stream,  considered  by  Ortloff 

2 

and  Ives,  can  be  described  mathematically  by  the  partial  differential 
equation, 


El  J  +  (M  +  m)  +  MU2  ix  +  2MU 

3x4  at^  ax*"  °'aX 


d_ 

ax 


_L  M  u2  (L-x)  & 
2  D  1  '  ax 


+  1  CN  S  U  H  +  U  S  -  0. 


(2.1) 


where 

El  =  bending  rigidity, 

M  =  lateral  virtual  mass  of  fluid  per  unit  length  of  wire  accelerated  by 
the  accelerating  wire, 
m  =  mass  of  the  wire  per  unit  length, 

U  =  velocity  of  the  free  stream, 

C-j.  =  drag  coefficient  due  to  pressure  acting  on  the  wire  surface, 

D  =  wire  diameter, 

L  =  total  wire  length,  and 

CN  =  drag  coefficient  due  to  shear  forces  acting  on  the  wire  surface. 


3 


TR  6343A 


The  special  case  of  a  zero  bending  rigidity  (or  an  infinitely  flexible) 
cylinder  is  realized  by  setting  El  =  0.  To  express  the  above  equation  in 
dimensionless  terms,  use 


r  =  (t/L)U  , 
e  =  L/D  , 


B  =  M/(M+m)  , 
n  =  y/L  . 


4  =  x/L  , 


The  above  equation  becomes 


2  ,2 

3  n  +  3  n 

1  ~~2  6 
3r  3| 


r  1  c  n  >n  +  an  <CT  +  CN> 
1  -  7  CT  6  (1  -6)J  Ti  - 2 


e8 


c«*6  I?-0- 


(2.2) 


The  associated  initial  boundary  conditions  are  described  by 


n  *  0 

|n |  is  finite, 
n  =  n1  (4) 


4=  0  (fixed  end  condition); 

4=1  (bounded  free  end  deflection); 
r=  0  (prescribed  initial  deflection); 
and 


(2.3) a 

( 2 . 3 )  b 

(2 .3) c 


an 

3  r  T=0 


0, 


0  (zero  initial  velocity). 


(2.3)d 


Ortloff  and  Ives  solved  the  problem  posed  by  equation  (2.2)  using 
conditions  described  in  equation  (2.3).  Their  solution  is  expressed  in  terms 
of  Bessel  functions. 
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The  initial  boundary  value  problem,  equation  (2.3),  for  the  partial 
differential  equation  (2.2)  is  said  to  be  well  posed  in  the  sense  of 
Hadamard'*  if  and  only  if  its  solution  exists,  is  unique,  and  depends 
continuously  on  the  data  assigned.  After  the  problem  is  formulated  using 
finite  difference  and  ordinary  differential  equations,  it  will  be  seen  that 
the  problem  is  well  posed.  We  will  seek  a  unique  solution  by  means  of  the 
numerical  techniques  presented  in  the  next  section.  When  the  boundary 
conditions  become  uncertain,  there  is  not  enough  information  available  to 
solve  equation  (2.2);  we  term  this  problem  ill  posed.  However,  a  derived 
boundary  condition  is  developed,  which  is  shown  to  be  adequate  for  our  problem. 

3.  THE  NUMERICAL  TREATMENT 

In  search  for  a  general  purpose,  accurate  solution  to  the  well  posed 
problem  (2.2),  subject  to  conditions  described  in  equation  (2.3),  the  method 
of  attack  is  to  discretize  ,  n^r  ,  and  by  central  and  backward  finite 
differences  and  then  to  transform  equation  (2.2)  into  a  system  of  second  order 
ordinary  differential  equations  (known  as  the  method  of  lines  ).  We 
discovered  that  Generalized  Adams  Bashforth  (GAB)  methods^’®  can  be  used  to 
solve  this  system  efficiently. 

Expressing  equation  (2.2)  in  short  form  and  writing  u  as  n  gives 

urf  +  a({)  +  buj  +  2s(u^)r  +  cur  -  0,  (3.1) 
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where 


am  =  B  [1  -|CTe(i-f)]  , 
b  =2^CT  +  CN)e®»  and 

c  =i8CTe. 

3.1  FINITE  DIFFERENCE  DISCRETIZATION 

Applying  the  second  order  central  and  backward  finite  difference 
discretization  to  equation  (3.1)  in  the  t  direction,  we  obtain 


u  . ,  -  2u  +■  u  ,  , 

/  \  .  / v  \  m+l  m  m-l  .  b  /  \ 

(u~>-  *  a(i->  - j? - *  F  lum  -  Vl> 


m 'tt  "nr 


*  ^  <“«  -  Vl>r  +  c<um>,  *  »• 


where  h  =  for  index  m  =  1,  2,  ... 

A  simplification  of  equation  (3.2)  gives 


(um)rr  +  Or  +  C)(um}r  *  TT  <  VlJT  +  ^  Vl 


(3.2) 


*b  2a(tm)1u  /b 

.h"  h2  JUm'\h  jum-l 


(3.3) 


6 


TR  6343A 


Equation  (3.3)  is  a  difference  equation,  representing  a  system  of  second 
order  ordinary  differential  equations  and  is  an  approximate  equation  to 
equation  (2.2). 


3.1.1  Consistency 

Before  we  apply  the  GAB  method  to  equation  (3.3),  let  us  examine  the 
consistency  of  our  finite  difference  operator  [u;h].  First,  expressing 
equation  (2.2)  in  true  operator  form,  we  obtain 


JC  [U] 


a(4)  b^ 

ar£  aC 


+  26 


3£3T  +  C  3T  ^ 


(3.4) 


Next,  expanding  u^  and  u^  in  powers  of  h  and  keeping  the  first  two 
principal  terms,  we  obtain 


Vl  =  um  +  h  (“m*'  +  T  (um)H  +  *  *  -  and 


Um  1  =  Um  “  h(Um)  ‘  +  (Uj" 

m-i  m  m  2  m 


Therefore, 


u  . .  -  2u  +  u  , 
m+1  m  m-1 


■  '“m’a 


I? 


h2{u»W 


(3.5) 
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and 


um  “  ur,i-l  a  (u„).  -  h2(u„)_ 
- R -  m{  m  {( 


(3.6) 


Substituting  the  power  expansions  of  u^  and  um_j  into  equation  (3.2) 
and  using  equations  (3.5)  and  (3.6),  we  find  that  equation  (3.2)  is  then 
expressed  in  a  difference  operator  form. 


£h[u;h]  = 


/  2 

2  .2  4  ' 

NVaU) 

a  .ha 
7  l?7 

\ar 

at  at 

+b 


a 

at 


ha2 


^7F 


a 

3 t 


h2  a2 


at 


+  cfelu  =  ° 


(3.7) 


* [u]  -  [u;h] 


.2  ,4  h  >2  h2  a2 

=  { -a( f ^  +  ^  b  — ^  +  26  j - 2 

“  ar  2  ar  i  at' 


It  is  seen  that 


lim 


UM  -  £h  MH])  >  °’ 


Therefore,  the  difference  operator  is  consistent  with  the  true  operator  in  the 
sense  of  Keller7.  Thus,  the  consistency  requirement  is  established. 
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Now  that  we  have  a  consistent  difference  operator,  we  seek  stable 
numerical  ordinary  differential  equation  methods  to  solve  equation  (3.2). 

3.2  ORDINARY  DIFFERENTIAL  EQUATION  SOLUTION 

To  seek  the  solution  to  equation  (3.2),  refer  to  equation  (3.3). 


Write  3u 

_ m  =  w 

3  T 


Equation  (3.8)  is  a  set  of  equations  that  represent  equation  (3.3)  as  a 
system  of  first  order  ordinary  differential  equations.  For  illustration, 
using  m  =  1,  2,  we  can  obtain 
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Uq  =  top  boundary  fixed  end  condition:  Ug  *  0,  £  ■  0; 

Wg  =  initial  condition:  uT  =  0,  r =  0;  and 
u-j  =  bottom  boundary  bounded  free-end  deflextion  condition: 
|u3|  is  finite  at  |=  1. 


The  matrix  elements,  A. of  matrix  A  can  be  determined  by  the  following 

*  J 

setups  in  which  we  define  A.  .  .  =  0  if  i-j  <  0  for  j  =  1,  2,  3. 

I  »  « 

When  index  i  is  odd,  A.  ^  =  1.  When  index  i  is  even, 

Ai , i-3  =  F"  ^  ' 
h 

rti , i-2  "  h  ’ 

a  2a(t  )  b 

Hi, i-1  "  h2  ‘  h  ’ 


where  a(|)  is  evaluated  at  a(SJ,  1=1 

l  2' 


Now,  the  problem  is  to  select  an  effective  numerical  ordinary  differential 


equation  method  to  solve  equation  (3.8).  A  close  examination  suggests  that 
the  Generalized  Adams-Bashforth  (GAB)  method  offers  an  efficient  solution.  In 
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the  present  application,  because  a  low  order  GAB  method  can  do  the  job,  high 
order  GAB  methods  are  not  necessary;  hence,  first  order  GAB  methods  were 
developed  into  computer  programs  in  FORTRAN  language.  However,  the  program 
package  is  flexible  so  that  high  order  GAB  methods  can  be  incorporated  when 
required. 

We  introduce  the  first  order  GAB 

n+1  Ah  n  .  .  .  ,BU . 

“  *  e  “  *  Ml,0<Ah>  \  (3.1 

to  solve  equation  (3.10),  where 

*l,0(Ah)  =  -(Ah)_1  (I  *  eAh)*  (3.] 

The  theory  with  respect  to  consistency,  stability,  and  convergence  has 
been  very  well  developed  for  Nonlinear  Multistep  (NLMS)  methods.®  The  GAB 
method  is  a  member  of  the  NLMS  family.  We  summarize  the  theory  below. 

NLMS  methods  take  the  expression. 


Jo  VAh<k~1)  u-*-  ■  h 


n*1  *  h  1^0  *ki  (Ah)  Vi 


3.2.1  Stability 

The  characteristic  polynomial  of  NLMS  is  defined  by 


(x,0  -  exAh  .iQ  a^1. 


(3.13) 


(3.14) 
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Using  the  GAB  method,  the  selection  of  oi  is  such  that 

ak  =  1,  <*k =  -1.  We  see  that  the  root  of  p(x,»)  has  unity  and  is 

simple;  therefore,  method  (3.11)  is  stable. 

3.2.2  Consistency 

The  GAB  method,  equation  (3.14)  satisfies  the  consistency 

condition 


lim 

h*0 


aie 


Ah(k-i) 

un+i 


h1J>  *ki(Ah)V( 


=  0 


(3.15) 


for  k  =  1,  «k  =  1  and  <*k_^  =  -1.  Therefore,  GAB  method  is  consistent. 

3.2.3  Convergence 

According  to  the  convergence  theorem  of  NLMS  methods,  "A  stable  and 
consistent  NLMS  method  is  convergent."  Therefore  the  GAB  method  applied  to 
problem  (3.10)  is  a  convergent  method. 

4.  BOUNDARY  CONDITIONS 


In  real  applications,  at  1,  the  bounded  free  end  deflexion  boundary 
condition  is  expected  to  be  such  that  n(l,»")  is  finite.  However,  the 
appropriate  function  n(l,T)  to  be  used  for  the  boundary  condition  appears 
uncertain  in  reality.  This  lack  of  information  defines  problem  (2.2)  as  an 
ill  posed  problem.  For  general  partial  differential  equations,  it  is  always 
difficult  to  formulate  correct  conditions  leading  to  a  well  posed  problem. 
Problems  may  look  reasonable,  yet  cannot  be  solved.  It  is  hoped  that  the 
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bounded  free  end  deflexion  boundary  condition  may  be  obtained  through 
experimentation,  but  the  exact  mathematical  expression  for  n(l,0  must  still 
be  worked  out.  We  will  attempt  to  change  the  ill  posed  problem  to  a  well 
posed  one  so  that  a  solution  exists  and  can  be  solved  by  the  numerical 
techniques  we  ha  *  developed. 

In  the  theory  of  second  order  partial  differential  equations  there  exists 
a  class  of  well  posed  problems,  such  as  the  Cauchy  problem  for  wave  equations, 
the  Dirichlet  condition  for  Laplace  equations,  and  the  mixed  initial  boundary 
value  problem  for  heat  equations.  Our  first  step  is  to  examine  the  most 
general  boundary  conditions.  Let  uN  denote  the  normal  derivative.  The 
first  boundary  value  problem  of  the  Dirichlet  type  indicates 

u  «  f  (4.1) 

on  the  boundary.  The  second  boundary  value  problem  of  the  Neumann  type 
indicates 


uN  -  f  (4.2) 

on  the  boundary.  The  third  boundary  value  problem  of  the  mixed  type  indicates 

uN  +  au  -  f  (4.3) 

14 
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on  the  boundary.  Note  that  the  third  boundary  value  problem  is  well  posed 
only  for  the  restricted  choice  of  a.  We  will  assume  that  the  free  end 
deflexion  boundary  conditior  takes  the  expression 

AU|y  +  au  =  f.  (4.4) 


When 


x  =  0,  a  =  1,  (4.4)  reduces  to  (4.1); 

A  =  1,  a  =  0,  (4.4)  reduces  to  (4.2);  and 

a  =  1,  a  arbitrary,  (4.4)  reduces  to  (4.3). 

In  our  application,  as  given  by  the  numerical  example  in  the  next  section, 

n(4,T)N  =  n(£,r)  ,  a  =  0,  a  =  1  gives 

n(|,T)  =  f(|,  )  and  Jf(l,  r  )|  is  finite.  This  gives  Ortloff's  and  Ives' 
bounded  free  encT  deflexion  boundary  condition. 

The  procedure  to  be  followed  here  for  determining  a  free  end  boundary 
condition  is  to  derive  an  approximate  boundary  condition  and  then  to  use  that 
boundary  condition  to  compare  the  solution  with  a  direct  application  of  the 
Ortloff  and  Ives  solution.  We  develop  a  form  of  the  boundary  condition  for 
the  second  order  partial  differential  equation  by  following  the  approach  used 
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by  Paidoussis  for  his  fourth  order  partial  differential  equation;  that  is,  by 
integrating  the  transverse  momentum  equation  over  a  short  tapered  end  which  is 
attached  to  the  free  end  in  order  to  generate  the  required  boundary 
condition.  Paidoussis  assumed  that  the  cross  sectional  area  tapers  smoothly 
from  S  to  zero  in  a  distance  (1)  sufficiently  short  that  the  forces  acting  on 
the  tapered  end  can  be  lumped  and  considered  in  appropriate  boundary 
conditions.  For  our  present  problem  the  boundary  condition  is  obtained  from 


L-I  (4.5) 


where  the  parts  of  the  equation  express  rate  of  change  of  fluid  momentum, 
hydrodynamic  forces,  and  cylinder  inertia,  respectively,  and  where  f  is  a 
factor  introduced  by  Paidoussis  to  account  for  the  intractable  flow  conditions 
at  the  free  end  and  v  is  the  transverse  velocity  of  the  fluid  relative  to  the 
cylinder.  Therefore, 


V  .  iZ  + 

at  ^ix’ 


(4.6) 
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F„  -  \  (5)  UCNV-  I4'7* 

T(x)  =  T(L)  +  \  (5)  U2Ct  (L-x);  (4.8) 

and  T(L)  is  a  consequence  of  form  drag  at  the  free  end. 

An  important  assumption  necessary  to  perform  the  integration  is  that  the 
length  of  the  tapered  section  (i)  is  small  enough  that  the  lateral  velocity 
(V)  may  be  considered  constant  over  1.  We  find 

«  f  ($  * "  4r)  * fuM  (»  * u  ») 

*  ($♦«$  *1  ©  "V® 

-T(L)«  4*71^t0lfZ)=°  (4.9) 

ax  at 

for  x  *  L,  all  t. 

After  nondimension  of  this  equation  as  before  and  neglecting  terms  of 
order  (|2)  and  t_  ,  we  have 

[f  *  £  (4)  £  *  [f  *  r  (4)  '  r  (4)]  It  ■ 0  t’-w) 

for  {  =  1»  al 1  T 
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On  physical  grounds  it  is  reasonable  to  neglect  Cy  relative  to 

g 

Cn/4,  making  the  final  boundary  condition 


an  + 


an 

ar 


=  0  for  Z 


1,  all  t, 


(4.11) 


which  amounts  to  a  "radiation  condition";  that  is,  no  reflected  energy 
exists.  In  the  following  sections,  we  refer  to  boundary  condition  (4.11)  as 
Kennedy  boundary  condition. 


5.  A  NUMERICAL  EXAMPLE 


The  test  example  is  obtained  through  linearization  of  a  fourth  order 
nonlinear  cable  equation,10  which  is  given  by 


(5.1) 


where  m,  M,  U  take  the  same  definitions  as  given  in  section  2.  T,  CN  are 
defined  as 


T 


(L  -  X), 


(5.2) 


T1  »  £T  w  D  U2  L, 


(5.3) 


18 


TR  6343  A 


111- 

CT  =  2  it  q2  lT,  and 


r  -  I  H  c  U 
LN  2  D  Lr* 


(5.4) 


(5.5) 


where  CN  and  Cy  satisfy  the  definitions  given  in  section  2. 


Assuming  a_  and  a_  commute,  expanding  equation  (5.1)  gives 
ax  at 

Performing  and  using  definitions  (5.2)  through  (5.6),  we  get 


P  -  7  CT  ^  (L-X)l  ‘“^(CT’CN) 

3 1  3X  J 


(5.6) 


.  on  3  V  +  1  r  i*  iZ  n 
+  2U  axat  7  CN  U  at  s 

(5.7) 


Equation  (5.7)  is  the  same  as  Ortloff’s  equation  (2.1)  before 
o 

nondimensionalization. 

Select 

D  -  r, 

2 

N  -  m  -  0.00273 
U  «  15  ft/sec 
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L  =  2000  ft 

Cj  =  1.8 

CN  =  1.1259. 


Then, 

e  =  48,000 


The  solution  to  equation  (5.7)  is  expressed  by 

y(t,  x)  =  e1wt  J7(x), 

where  7  is  approximately  21.89  and  J^x)  initial  values  are  calculated  using  a 
UNIVAC  1108  Bessel  function  subroutine.^ 

The  fixed  end  boundary  condition  initially  is  zero.  The  free  end  boundary 
condition  uses  n(l,*r)  =  e1WTJ7(l). 

This  problem  was  tested  again  using  Kennedy's  free  end  boundary 
condition.  Results  are  surprisingly  in  agreement  with  the  known  solution. 

The  test  results  seem  to  show  that  the  Kennedy  free  end  boundary  condition  is 
adequate.  Results  are  presented  in  graphic  form.  Two  sets  of  graphs  are 
given:  one  displays  M£,t)I  versus  t,  the  other  displays  the  real  {n(|,r)f 
versus  t.  Both  plots  are  constructed  at  §  *  0.2,  0.8. 
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{  =  0  2 
.5 


.5 _ 

6  10  20  30  40  50  60 

TIME  *  102 


Figure  1:  Solution  magnitude  vs  time 


£  =  0  8 
.5 

I  A 
1^,.  -I 

.5i - - - 

0  10  20  30  40  50  60 

TIME  *  102 


Figure  3:  Solution  magnitude  vs  time 


{  =0.2 


0  10  20  30  40  50  60 

TIME  *102 


Figure  2:  Real  part  solution  vs  time 


{  =  0.8 


w  0  10  20  30  40  50  60 

TIME  *102 


Figure  4:  Real  part  solution  vs  time 


21 


TR  6343A 


6.  CONCLUSIONS 

A  numerical  solution  to  the  dynamic  motion  of  a  zero  bending  rigidity 
cylinder  in  a  viscous  stream  has  been  introduced.  The  numerical  procedures 
developed  to  obtain  the  solution  are  theoretical ly  convergent  and 
computationally  accurate. 

For  given  appropriate  boundary  conditions  and  accurate  initial  values, 
this  model  will  produce  an  accurate  unique  solution.  For  uncertain  boundary 
conditions,  this  model  can  be  used  as  a  tool  to  study  the  boundary  effects  and 
possible  to  construct  the  ad  hoc  boundary  conditions. 


22 


TR  6343A 


REFERENCES 

1.  M.  P.  Paidoussis,  "Dynamics  of  Flexible  Slender  Cylinders  in  Axial  Flow," 
Journal  of  Fluid  Mechanics,  vol.  26,  part  4,  1966,  pp.  717-751 

2.  C.  R.  Ortloff  and  J.  Ives,  "On  the  Dynamic  Motion  of  a  Thin  Flexible 
Cylinder  in  a  Viscous  Stream,"  Journal  of  Fluid  Mechanics,  vol.  38,  part 
4,  1969,  pp.  713-720. 

3.  P.  R.  Garabedian,  Partial  Differential  Equations,  John  Wiley  3  Sons, 

Inc.,  NY,  1964. 

4.  W.  F.  Ames,  Numerical  Methods  for  Partial  Differential  Equations, 

Academic  Press,  NY,  1977. 

5.  D.  Lee  and  S.  Preiser,  "A  Class  of  Nonlinear  Multistep  A-Stable  Numerical 
Methods  for  Solving  Stiff  Differential  Equations,"  Journal  of  Computers 
and  Mathematics  With  Applications,  vol.  4,  no.  1,  1978,  pp.  43-51. 

6.  D.  Lee,  Generalized  Adams  Methods,  NUSC  Technical  Report  No.  6011,  Naval 
Underwater  Systems  Center,  New  London,  CT,  1979. 

7.  H.  B.  Keller,  Numerical  Methods  for  Two-Point  Boundary-Value  Problem, 
Blaisdell  Publishing  Company,  Waltham,  Massachusetts,  1968. 

8.  D.  Lee,  Nonlinear  Multistep  Methods  for  Solving  Initial  Value  Problems  in 
Ordinary  Differential  Equations,  NUSC  Technical  Report  No.  4774,  Naval 
Underwater  Systems  Center,  New  London,  CT,  24  May  1974. 

9.  S.  F.  Hoerner,  "Fluid-Dynamic  Drag,"  published  by  the  author,  1965. 

10.  J.  Patel,  personal  communication,  1979. 


23/24 
Reverse  Blank 


TR  6343A 


APPENDIX 

COMPUTER  PROGRAMS  STRUCTURE  AND  COMPUTER  LISTING 

COMPUTER  PROGRAM  STRUCTURE 


MAIN 
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ACRONYMS 

MAIN  main  program  which  controls  the  setup  of  inputs  and  the  preparation  of 
outputs 

START  supplies  the  initial  values 

DIFEQ  controls  the  present  T-step  and  calls  for  NLMS(GAB)  method 

NLMS  1st  order  Generalized  Adams-Bashforth  method 

INVER  calls  for  complex  matrix  inversion 

C6JR  complex  matrix  inversion  using  Gauss-Jorden  reduction 

GFN  calculates  the  g-vector 

BC  fixed  end  and  free  end  boundary  conditions 

PADE  a  rational  function  approximation  for  matrix  exponentials 

INVERT  calls  for  double  precision  matrix  inversion 

DGJR  double  precision  matrix  inversion  using  Gauss-Jorden  reduction 

The  user  needs  to  deal  with  MAIN,  START,  GFN,  and  BC.  The  user  need  not 
be  concerned  with  DIFEQ,  NLMS,  INVER,  INVERT,  CGOR,  and  DGJR. 
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COMPUTER  PROGRAM  LISTING 


1 


MAIN 


COMMON  PEPS  *  PCT . PCN , PBETA , EQ , PA , PB  , PC , PD  , DZ , IBND 
COMMON  A( 18, 18) ,T(3) 

D T  MENS  TON  Y ( 2  *  1 8 ) , YZERO (18)* YNEW (18)*  EXACT (18) 

COMPLEX  A  *  Y  *  YZERO , YNEW  *  XX  *  OMGA *  SAVE  *  EXACT 
*****  THIS  PACKAGE  SOLVED  A  2ND  ORDER  P.D.E.  BY  THE  METHOD 


OP  LINES  AN 


***** 

***** 

***** 

***** 

***** 

***** 

***** 

***** 


***** 

***** 

***** 

)l(  )|(  )|[  )|( 


DEFINITIONS  * 


G E N E R A L I Z E D  ADAMS- BASH E 0 R T H  M E T W 0 D S 

reference;  ortloff  and  ives 

INPUT  PARAMETERS  HAVE  THE  FOLLOWING 
N  ~  NUMBER  OF  2ND  ORDER  ODE 

TMAX  -  MAXIMUM  TAO 

TINT  =  EUERY  TAO  INTERNAL  TO  BE  PRINTED  OUT 
PXI  =  AT  THIS  XI «  THE  OUTPUT  IS  REQUESTED 

FQ  •-  FREQUENCY  ?  NOND I  MENS  TONAL  OMEGA 

IBND  =  BOUNDARY  CONDITION  INDICATOR 

=  1  BUILT-IN  KENNEDY  BOUNDARY  CONDITION 
■  2  USER-SUPPLIED  BOUNDARY  CONDITION 
=  EPSILON 


PEPS 

PCT 

PCN 

PBETA 

H 

READ 


=  C  SUB  T 
=  C  SUB  N 
=  BETA 
=  TAO  STEP 
INPUTS  HERE 


SIZE 


N, TMAX, TINT, PXI.FG, IBND 
PEPS , PCT , PCN , PBETA , H 


READ (5,*) 

READ (5,*) 

N-2#(N-1 ) 

DZ-l ,0/( FLOAT (N/2+1 ) ) 

PB=0 . 5*PBETA*PCT*PEPS 
PA=PBET  A*  < 1 . 0-PB ) 

PD=0 . 5*PCN*PEPS*PBETA 
PC=PB+PD 

*****  TO  SET-UP  MATRIX  A 
DO  20  1=1, N 
DO  20  J=1,N 
20  A ( I , J) =CMPLX (0.0,0 
BH=PC/DZ 
TBH=2  *  #PBETA/DZ 
DG=- ( TBH+PB ) 

DO  28  1=1, N 
IF ( I  . GE *  4)  GO  TO 
IF (I  . EQ .  2)  GO  TO 
A ( I ,  I  + 1 ) =CMPLX ( 1 . 0 , 0 ♦ 0 ) 

GO  TO  20 

25  J=I/2 

X=PBETA* ( 1 . - . 5*PCT*PEPS* ( 1 . - J*DZ > ) 
A(I,I-1)=CMPLX(2.*X/(DZ*DZ>-BH,0,0> 
A  ( I ,  I )  =CMPLX ( DG , 0 . 0 ) 

A (1,1+1 )=CMPLX(-X/(DZ*DZ) ,0.0) 

GO  TO  28 

26  IF (MOD (I, 2)  .EG.  0)  GO  TO  27 
A(I,I+1 )=CMPLX(1 .0,0.0) 

GO  TO  28 


.0) 
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-PBETA* 


K' 


i  .-jkDz; 


*  ?2i  0 . 0 


c  A  v  E  =  Y  Z  E  ft  0  \  M  -  i  > 

10  CONTINUE 
TX-TX+TINT 

[:ALL  D I F E  0  ••  H  ?  m  ?  T  X  »  Y  ?  Y 
C  *****  RESULTS  IN  '<  H  E  W  (  I  > 

•1'  *****  Z  CONTAINS  PRESENT 
C  t  *  4  f  .*  U  SER  J  5  E  3  A  3  0 v  E  IN  I" 
SAVE  =  YNEW (  N~  1  ) 


R0>Z>YNEW»3A V E ) 


3 


2  F  T  7 1 P  T  -i-  i 
IFCFT  , L£ 

irirtoiK  r-  - 

SO  TO  100 
CONTINUE 


4)  SO  TO  3 

.EG.  0 >  00  T  0  3 


WRITE (5fl)  H»Z 

1  FORMAT  (  / 1  0  X  >  '  H  -  '  ?  E 1 5 . 3  *  5  X  .*  '  T  =  '  ,  E  1  5 . 8 

2  FORMAT  <  3X  f  F3 ♦ 2  r 1  OX  f  E15  *$  » lOX  ?  E15 »  S  > 

0  *****  PRINT  OUT  EXACT  SOLUTION 

C  *****  THIS  F  3RTI0N  13  FOR  TEST  EX  A  M F  L E  0  N  L Y 
CALL  START (N»T  »  E  X  A  C  T  > 

DO  100  I  =  1 f  N  f  2 
NN  = ( I + 1 >/2 


D  L  Z  =  N  N  *  D  Z 

WRITE <5, 2) 

DLZ  t 

Y  NEW  (  I  ) 

WRITE(5f3) 

EXACT 

I  ) 

3  FORMAT  (  2 1 X 

r  E 1  5  .  S 

>10X»E15.3/> 

100  CONTINUE 

IF  ( TX  .GE. 

T  M  A  X  ) 

STOP 

GO  TO  10 
END 
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NLMS 

SUBROUTINE  NLMS ( H r v  * N ♦ YN t IS* SAVE > 

COMMON  PEPS  *  PCI  *  PCN y  PBETA *  FQ *  PA t PB , PC  *  PD * DZ » I BND 
PARAMETER  KM=l8 
COMMON  A  ( KM ?  KM >  *  T  (  3 ) 

DIMENSION  AH ( KM  * KM ) *  EAH ( KM  *  KM ) *  G  <  KM ) 

D I MENS I ON  P I ( KM  *  KM ) *  UN I T < KM  *  KM > 

D I MENS I ON  Y  ( 2  * KM ) y YN ( KM ) 

D I M E N S 1 0 N  F E <  KM  *  KM ) ?  A1 <  KM  * K M ) 

COMPLEX  A  v  AH y EAH * G  r PI * PH  r  UN I T  y Y  r  YN  *  SAVE 
COMPLEX  FE  y  A1 
DATA  IND/O/ 

IF ( IND  *  G  T «  0)  GO  TO  14 
DO  2  1*1  yN 

DO  2  J“ 1 y  N 

2  AH<IyJ)=H#A(I»J) 

TF(N-l)  7*8*7 

8  A:l.  <  1.  *  1  )-l  ,/AH(l  y  1  ) 

GO  TO  9 

7  CALL  INVER (AHy Ny Al ) 

9  IFCN— t )  1 0 . 1  •:  *  I  0 

1. 1  EAH  <  t  y  1. )  --CEXP  (  AH  ( 1.  y  1 )  ) 

GO  TO  1 4 

3.0  CALL  PA  DE  ( A  *  l-l  *  EAH  *  N  > 

IND  =  IND  +  :l 

14  IF (IS  ,Gi.  I)  GO  TO  100 
DO  •  I - 1  *  N 

DO  A  J*t *  N 

P 1  ( I  y  J 1  ~CMP!.. X  ( 0 , 0  y  0 , 0  > 

UNI T  ( I  *  J  > --CMPI.  X  ( 0 ,0*0, 0 ) 

6  CONTINUE 

1  UNIT ( I y I ) -CMP LX ( 1 . 0*0, 0  > 

100  CONTINUE 

CX  ^  ^  \L«  ^  O/  ^  ^  «J/  ^  ^  vb  U/  ^  vLr  ^  ^  ^  U#  ^  ^  ^  J*  ^  Oy  ^  ^  .  .X»  ■!.  ^  J,  ^  ^  ^  ib  ^  .  i  .  > L,  ^  ^  .b  ^  ^  'b  ^  ^  ^  ^  ^  »b  >b  ^  ■  b  ^  ^ 

*  ^  ^  ^  *T  ^  *  -T  ^  ™  ^  ^  yp  •*  *  ^  ^  *•  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ®  ^  ^  *•  ^  ^  ^  ^  ^  ^  ®  ^  ^  ^ 

C  *  NONLINEAR  MULTISTEP  STARTS  HERE,  * 

C  *  BEGINNING  SECTION  DOES  INITIALIZATION  * 

0  j**#*  ********  t#%**%)l()l()l<xx*x*XX*)l(yxx********#x*#*)K 

DO  I 32  1*1 fN 
1  32  YN  ( I )  -CMPI.  X  (  0 ,0*0.0) 

IF  ( IS  ,  GT ,  1 )  GO  TO  1.3! 

DO  103  1*1.  *N 
DO  103  J*1.  *N 

103  PI  <  T  *  J)*-EAH(  T  *...1)  fUNTT‘  I  yj) 

C  *$**  ********  #***  ******%>i(**####**XX#**##****x 

-  *  1ST  ORDER  GAB  * 

C  *  DO  LOOP  105  CAI.C, HATES  PHI  (1*0) 

C  *  LOOP  108  OR  11.0  COMPUTES  FINAL  Y(N+1) 

***  #*  *  *  *******  ***************  *  **  **  *****  *****  *  *****  *  ******  **  ***  ***  * 
on  i  oa  i; - 1 ,  n 
DO  106  .  J— 1  * N 
PH*  CMP"  I  X(  0,0,0.  O' 

DO  105  K* I *N 

I  05  AH--PH-A1  ( I  * K > *P1  <K*.D 
FE ( T  *  I ) -PH 
1 08  CONTINUE 

*  Vi  CAL  I  OF*  NIG*  H  y  N  *  Y  *  1  *  I  *  A  *  GAU.*-  > 
i  n  i  *i  0*8  J  *  i  *  N 

no  1 08  ,i*i  »n 

i 08  YN ( T )*YN f I ' +FAH (I*  n  * Y ( 1  ,  1 1 +H*FF O .  i  >  *G (  J  ) 

rftupn 
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DO  9  1=1, N 

DO  10  J=1 , N 
C  <  I ,  J ) =AA  <  I » J  >  #H/2 . 0 
PPdf  J)=O.DO 
10  CONTINUE 

C(I,I)=C(I, I)+1.D0 
9  CONTINUE 
DO  12  1=1, N 
DO  13  J=1 ,  N 
DO  14  K=1,N 

P  P  (  I ,  J )  =  P  P  ( I ,  J )  +  B  ( I ,  K  )  *  C  ( K ,  J ) 

14  CONTINUE 
13  CONTINUE 
12  CONTINUE 

IF < M  .EC).  0)  GO  TO  40 

^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  j|^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^  ^ 

C  *  NORM (AH).GT.(.l),  EXP(A>  -EXP < A/2*#M ) ** ( 2**M )  * 

DO  24  1=1, N 

DO  25  J=1 , N 
B ( I , J ) =0 . DO 
25  CONTINUE 

24  CONTINUE 
DO  36  K=1 , M 
DO  27  1=1, N 
DO  28  J=1 , N 
DO  29  L=1,N 

B ( I ,  J )  =B <  I ,  J )  +PP (1,1..) *PP (L,  J) 

29  CONTINUE 
28  CONTINUE 
27  CONTINUE 
DO  31  1=1, N 
DO  32  J=1,N 
BB=B ( I , J ) 

P(I, J)=CMPLX(BB,0.0) 

B  < I , J  >  =0 . DO 
32  CONTINUE 
31  CONTINUE 
36  CONTINUE 
H=HAVE 
RETURN 
20  H=H/2 . 0 
M=M+1 

DO  54  1=1, N 
DO  55  J=1 , N 
PP ( I , J ) =0 . DO 
55  CONTINUE 
54  CONTINUE 
GO  TO  30 
40  H=HAUE 
RETURN 
END 
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START 


SUBROUTINE  START ( N t T p YZERO ) 

COMMON  PEPS  p PCT » PCN  t PBETA »  FQ  p PA  * PB  t PC p PD  t DZ 1 1 BND 
COMPLEX  YZERO (1) 

DIMENSION  T < 1 ) 

******  USER  SHALL  REVISE  THIS  PORTION  TO  INCORPORATE  HIS  INITIAL  VALUE 


♦  ♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 


YZERO ( I )  CONTAINS  THE  FUNCTION 

YZERO <1  +  1 )  CONTAINS  THE  DERIVATIVE 


DO  29  L.P=1pNp2 
DLZ=DZ*  <  LP+1 ) /2 

IF  <  LP . EG  *  1 )  YZERO  <  LP ) =CEXP  <  CMPLX  <  0 . 0  p FQ*T < 1 ) ) ) 

*  *CMPLX <  *  9954A739 p - . 73021506E-01 ) 
IF ( LP  *  EQ ♦ 3 )  YZERO  <  LP ) =CEXP<  CMPLX ( 0 . 0  p FQ*T ( 1  >  ) ) 

*  *CMPLX( .98A02435p-. 14547502) 

IF  <  LP . EQ . 5  >  YZERO  <  LP ) =CEXP <  CMPLX  <  0 . 0  p FQ#T < 1 ) ) ) 

*  #CMPLX< . 971 72089 p -♦ 21 705468 > 

IF  ( LF'  *  EQ » 7 )  YZERO ( LP ) =CEXP < CMPLX <  0 . 0  p FQ#T ( 1 ) ) ) 

*  *CMPLX( . 96262427 p -* 28745766 ) 
YZERO  <  LP+1 )  '--YZERO  <  LP  >  *CMPLX  <  O.OpFQ) 


29  CONTINUE 
RETURN 
END 


GFN 


S LI r  R 0 UTI N E  G F N  (  3  p  H  t  N  r  Y  ?  J 


A  i  S  A  V  E  ) 


C  *  *  *  *  *  THIS  HANDLES  THE  Q  VECTOR 
C  *  +  tc  +  ★  G  VECTOR  CONTAINS  BOUNDARY  INFO F; M A T 1 0 N 

COMMON  F  Er  5  p  F  CT  r  PCN  r  FBETA  p  FG  p  F'A  f  FB  p  F'C  p  F  D  p  DZ  p  I BND 
C  0  M  F  L  E  X  A  <  1 S  f  1 S )  ?Y(2f13)  »6(1S)  f  SURF  j  BOTT  f  XA  p  XE:  f  SAVE  r  WZ 
D I h E N  SIGN  Til) 

DATA  PI f ZERO f ONE/3. 1415926535 fO. Of  1 .0/ 


J  C‘  1  I  ~  1  f  N 

.i  O'  I )  =  C  H  F'  L  X  i  0 . 0  f  0  «  0 ) 

C  *  +  *  *  *  FIRST  ARGUMENT  Of  CALLS  FOR  FIXED  END  CONDITION 
CALL  3  C  <  0  p  N  f  H  f  T fY  fSAVEfWZpSURF) 
u  <  2  )  -  F  B  E  T  A  *  (  1 ,0-0.5*PCT*PEPS*<  1 . 0-DZ>  ) 

G  2  )  r  2  .  G*r  BETA  +  WZ  / DZ  +  -  F'C/DZ-G  t  2  >  /  DZ#*2  )  *5URF 
W  RITE  i  5  f  *  )  6  ■"  2  ) 

*****  FIRST  ARGUMENT  If  CALLS  FOR  FRSE-ENH  CONDITION 
CALI.  BCi  IfNfHfT  fY  fSAVEfWZfBOTT) 

X  =  -F  BETA * ( 1 . - . 5*PCT*FEPS* <1 ♦-< N/2) *DZ ) > / < BZ*DZ ) 

G  ( N )  --CMFL.X  <  X  f  0 . 0  )  *BOTT 
WRITE- 5f * -  Si N) 

RETURN 

HMD 
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INVERT 


SUBROUTINE  INVERT <A»N»ANS) 

****************************************************************** 

*  MATRIX  INVERSION  SUBROUTINE ,  CALLED  BY  FADE  OR  NLMS  * 

*  A  CONTAINS  THE  ORIGINAL  ELEMENTS  AND  REMAINS  UNALTERED  * 

*  ANS  CONTAINS  THE  A**<-1>  * 

*  THIS  SET-UP  IS  USING  UNIVAC  1108  MATHF’K  EXISTING  DOUBLE  * 

*  PRECISION  GAUSS-JORDAN  REDUCTION  * 

*  THIS  PROGRAM  IS  REPLACEABLE  BY  THE  USER  * 

****************************************************************** 
DOUBLE  PRECISION  A< 18r 18) >ANS< 18> 18) >V(2) 

DIMENSION  JC< 18 ) 

DATA  NR/18/»NC/18/ 
v<n=i.Do 
DO  1  1=1 fN 
DO  1  J=lfN 

1  ANS  ( I » J )  =A  <  I » J ) 

CALL  DG JR ( ANS  t NR » NC  t N » N » MDEX  r JC  t V ) 

IF < MDEX  .EQ.  1)  GO  TO  10 
RETURN 

10  WRITE  <  4 » 2 ) 

2  FORMAT (3X,22HMATRIX  INVERSION  ERROR) 

RETURN 

END 


INVER 


SUBROUTINE  INVER < A t N » ANS ) 

PARAMETER  NDIM=18 

COMPLEX  A(NDIM,NDIM) >ANS(NDIM.NDIM) fV<2) 
DIMENSION  JC(NDIM) 

DATA  NR/NDIM/rNC/NDIM/ 

V(1)=CMPLX(1. 0,0.0) 

DO  1  1  =  1  »N 
DO  1  J=1 t N 

1  ANS( I » J)=A< I » J) 

CALL  CG JR  <  ANS » NC . NR » N > N > MDEX  >  JC»  V ) 

IF < MDEX  .EQ.  1)  GO  TO  10 
RETURN 

10  WRITE  (4.2) 

2  FORMAT ( 3X» 1 1HMAT  INV  ERR) 

RETURN 

END 
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DGJR 


SUBROUT I NE  DGJR < A » NC » NR > N > MC > MDEX >  JC  >  U ) 
DIMENSION  JC<N)»U<2> 

DIMENSION  A<NR»NC) 

DOUBLE  PRECISION  A»U»X»DLOG 
IW=V<1> 

M=1 
S=1 . 

L=N+<MC-N>#<IW/4> 

KD=2-  MOD ( IW/2  r  2 ) 

IF <KD*EQ» 1 )  U  <  2 ) =0  * 

KI=2-  MOD < IW » 2 ) 

GO  TO  <5*20) j  K I 
5  DO  10  1  =  1  *N 

10  JC<I)=I 

20  DO  91  1  =  1  *  N 

GO  TO  <22*21 >  *KI 

21  M=I 

22  IF  <I.EQ,N>  GO  TO  60 

X=-l  ♦ 

DO  30  J=I  *  N 

IF  <  X ♦ GT ♦ ABS <A<J*I))>  GO  TO  30 
X=ABS  <A<,J*I)> 

K-J 

30  CONTINUE 

IF(N.EQ.I)  GO  TO  60 

U<1)=-U<1) 

GO  TO  < 35 » 40 ) r NI 
35  MU= JC ( I ) 

JC  < I ) =JC  <  K ) 

JC  <  K )  =M(J 

40  DO  50  J=M.«L 
X=A< I  *  J) 

A  <  I  *  J  )  =A  <  K  *  J ) 

50  A  <  K  *  J  > =X 

60  IF  ( ABS (A(IfI))« GT ♦ 0 * )  GO  TO  70 
IF < KD *EQ ♦ 1 )  U<1>=0* 

JC(1)=I-1 

RETURN 

70  GO  TO  <  71  *  72 ) * KD 

71  IF<  A< I  *  I > .LT ♦ 0. )  S=-S 

0<2)=U<2)  +  DLOG<ABS<A<  1 * 1 )) > 

72  X=A < 1 * 1 ) 

A<I*I)=1» 

DO  80  J=M * L 

A  < I  *  J ) =A  < I  *  J ) /X 
CALL  ERRTST  <  72  * MDEX ) 

IF<MDEX «  EQ» 1 )  GO  TO  150 

80  CONTINUE 
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DO  91  N=1»N 

IF  (K.EQ.I)  GO  TO  91 
X=A(Kj»  I  > 

A  <  N » I  >  =0 . 

DO  90  J  =M>L 
A  <K » J> “A  <K»  J  > -X#A ( I » J  > 
CALL  ERRTST <72»M0EX> 
IF<MDEX.EQ. 1 >  GO  TO  150 

90  CONTINUE 

91  CONTINUE 


GO 

TO 

(95 

» 140)  i 

>  KI 

DO 

130 

J- 

1 » N 

IF  < 

JC< 

J)  . 

EQ. 

J) 

GO 

TO 

130 

JJ= 

‘J+l 

DO 

100 

1  = 

JJ» 

N 

IF  <  JC ( I > . 

EQ. 

J) 

GO 

TO 

110 

100  CONTINUE 
110  JC<I)=JC(J) 

DO  120  K=1»N 
X=A ( K  f I ) 
A(KrI)=A(KfJ) 

120  A ( K»  J) =X 
130  CONTINUE 
140  JC  < 1 ) =N 

IF ( KD ♦ EQ ♦ 1 )  U<1)=S 
RETURN 

150  JC(1)=1-I 

IF < KD . EQ . 1 )  V<1)=S 

RETURN 

END 
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CGJR 


SUBROUTINE  CGJR <A»NC»NR,N,MC,IFL,JC»U) 
DIMENSION  JC<1) 

COMPLEX  CL„0G,U,XC,A<18,18> 

COMPLEX  7 

INTEGER#2  NERR 

NERR=72 

IFL=0 

IW=U 

U=<0. ,0. ) 

IBIT=0 

M«1 

L=N+ ( MC-N )  #  < IW/4 ) 

KD=2-M0D < IW/2 » 2  > 

KI=2~M0B < IW , 2 ) 

GO  TO  ( 5  j>  20 )  ,KI 
5  DO  10  1  =  1, N 

10  JC<I)=I 

20  DO  91  1=1, N 

GO  TO  <22,21 ) ,NI 

21  M=I 

22  IF  <I.EQ.N>  GO  TO  60 
X=- 1 . 

no  30  J=I,N 

ANORM=ABS  <  REAL  <A(J,I))) +ABS ( AIMAG  <A(J,I>) ) 

IF(X.GT *ANORM)  GO  TO  30 

X=ANORM 

K=J 

30  CONTINUE 

IF(K.EQ.I)  GO  TO  60 
IBIT=IBIT+1 
GO  TO  <35,40) ,KI 
35  MU= JC  < I ) 

JC  < I ) = JC  <  K ) 

JC  <  K ) =MU 

40  HO  50  J=M,L 

XC=A  < I , J ) 

A  < I , J ) =A  <  K , J ) 

50  A  <  K , J > =XC 

60  ANORM= ABS<REAL<A<  I , I ) ) ) +ABS < AIMAG < A < I , I ) )  ) 

IF  <  ANORM ♦ GT  *  0 )  GO  TO  70 

U=<0, ,0. ) 

JC<  J.  )  =  I-1 
RETURN 
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70  GO  TO  (71*72) >KD 

71  U=U+CLQG(A<I.I> > 

Z-CLQG ( A ( I »  I )  ) 

72  XC=A/IfI) 

A(I»I)*(l»iOt  ) 

DO  SO  J»M,L 
A(IfJ)=A(IrJ)  /XC, 

CALL  ERR TST<NERR»IFL> 
1F<IFL,EQ.1)  GO  TO  150 

80  CONTINUE 

DO  91  K=1 *N 

IF  ( K • EQ . I >  GO  TO  91 

XC-A<K* I ) 

A<K.I)=(0* .0.  ) 

DO  90  J  =Mf L 
A  <K»  J )=A<  K  .  J ) -XC#A  < I » J) 
CALL  ERRTST (NERR. IFL) 
IFIIFL.ECM)  GO  TO  150 

90  CONTINUE 

91  CONTINUE 

GO  TO  <95.140) .KI 


DO 

130  J 

=*1 » N 

IF 

< JC< J) 

.EQ. 

J) 

GO 

TO 

130 

JJ= 

=  J+1 

DO 

100  I 

-  J  J » 

N 

IF 

<  JC  < I  ) 

♦  EG. 

J) 

GO 

TO 

110 

100  CONTINUE 
110  JC<I)=JC<J> 

DO  120  K-l rN 
XC=A<K.I> 

A  <  K r I ) =A  <  K . J ) 

120  A <K.J)-XC 
130  CONTINUE 
140  JC  < 1 >  «N 

U=U+<0.  .3.14159265)*CMPI,X<FL0AT<M0D<IBIT,2)  )  .0,  ) 
RETURN 

150  JC<1)»1-I 

RETURN 
END 
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BC 


•j  ?  :\*  ?  n  '  :  ? 


7  *  i  *  7 


y  z  i“!  V  Z.  ?  AJ  ,V  ?  V  ; 


«■  7  *.  Y  Y  L'SEF  h:ic  AN  :.:r7I0N  TO  £-ir  F  L  HIS  BOUND ARY  03  Pit  I 
MnJ  I THI3  CASE  *  57 aTEMENT  PF7ER  5  5HC-U:_D  5 E  F,  EFLF 
CONE  I  7  I  0  N 

i-  G rt i“  J N  r‘ L F' :!  f  FCT  tP C N  ?  F  r  E  f.A  t  F Q  ,■  F'  A  »  F  S  7  F' C  7  F' D  »  £' Z  ?  I  r  N li 
0  I  h  E  N  3 1 0  N  T  ■,  :■ 

1  0  N "  L  E  X  Y  \  2  »  i  5  ?  a  7  \j  '< 

Ml*<  r  L  X i.'  E  l.«  r.i C :J  N  D  A R  1  G N Dili  C N 
’  •  E  X 3  H F  i_  X  .  j  .  0  - G  f  7  •  1  ■  ;■  '• 

■i  .<  =  I-  *  F’ L  X  i  0  .  v  >  F  0  )  !  ■; 

E  1  j 

•  • :  :  ...  >  r  l  !j  i; 

.1  7  3  :  •  2  7  3  ■  *  1 5MD 

;  1  m  1  *vk*  »  n.  -•  M  *  f  t  ‘  ♦  t- 

7  7  Z  t- 

■  *  •?  5  >  t  iii*M  * h  U * i *  7 

C  G  N  7 1  rJ  Z  E 

>'  *  <  H  k '1  (  1  7  N- 1  >  D Z r 3  A v E  )  /  <  i  ,  0 r  H / D Z  ) 

RETURN 

r*  n*  USER  SUPPLIED  BOUNDARY  CONDITION 
2  CONTINUE 

a  =CSX F  {  CMF  LX  (  0 . 0  »  F Q * 7  v  1  >  )  }  kCMPLX  (  .  9 2 ?  3 6  i  ? 6  ?  -  .  E 3 c 7 
X-ChPLA< i.OfO.O) 

RETURN 
2  CONTINUE 

I F  FEET  A  *  N  E  ,  1  <  0 )  30  TO  21 
X  ~  Y  ( 1  ?  N  - 1  > 

REIT  URN 

31  V:Z*1 .  +  •  PCT:*PEPSkPBETA*IiE>/(2,*<  1  .  -F'BETA  )  ) 

X  -  <  C  M  F  L  X  (  1  .  /  X  Y  Z  f  0 . 0  >  )*<GMF'LX<  1  .  +  X  Y  Z  1 0  •  0  >  #  Y  < 1 »N-1 > 

RETURN 

END 


0  7  ?  9  > 


-  1  iiilrj 
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SUBROUTINE  DIFEQ ( H * N » TMAX , Y » YZERO t ANORM , YN > SAME ) 

COMMON  PEPS t PCT , PCN t PBETA » FQ , PA »PB»PC» PD > DZ * IBND 
COMMON  A ( 18  » 18 ) *  T  <  3  ) 

COMPLEX  Y(2»18)kYN(18>» YOLD ( 2 » 18 ) * YZERO < 1 8 ) » A » SAVE 
DO  40  1=1 fN 
40  Y  < 1 » I ) =YZERO  < I ) 

IH=0 

TZERO=T ( 1 ) 

60  TEA=T ( 1 ) +H 

IF ( TEA . GT . TMAX )  H=TMAX~T ( 1 > 

IF  <  TEA . GT . TMAX )  IH=0 
IF  <  TEA . GT  *  TMAX )  GO  TO  60 
IH-IH+1 

IF < IH  .GE,  32767)  IH=2 

T ( 2 ) =T  < 1 ) +H 

IMP=2 

DO  62  J=1 r IMP 
DO  62  1=1. N 
62  YOLD ( J  *  I ) =Y  <  J  f I ) 

CALL  NLMS ( H , YOLD . N » YN » IH » SAME > 

59  DO  66  1  =  1  .N 
Y ( 2 » I >  =YN <  I ) 

66  YOLD (2*1) =YN ( I ) 

Xc***************************************************************** 

*  RESULTS  Y ( TEA )  IN  YN<I>  AND  Y(2fl)  * 

****************************************************************** 

ANORM=TEA 
DO  85  I=l« N 
Y ( 1 » I ) =Y  <  2 » I ) 

YZERO < I >  =Y< 1 » I ) 

85  CONTINUE 
T ( 1 ) =T ( 2 ) 

TZERO-T <1 ) 

IF ( ABS  <  TEA-TMAX ) . LE ♦ < . IE-5) )  RETURN 

GO  TO  60 

END 
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PADE 

SUBROUTINE  PADE <AjH»P»N> 

C  *****  A  RATIONAL  APPROXIMATION  OF  MATRIX  EXPONENTIALS 
C  *****  DOUBLE  PRECISION  IS  NEEDED  FOR  REQUIRED  ACCURACY 
PARAMETER  NM=18 
COMPLEX  A ( NM  t NM ) * P ( NM , NM ) 

DOUBLE  PRECISION  AA ( NM , NM ) » PP ( NM t NM ) » B ( NM t NM ) » C ( NM > NM ) > HAVE 
DOUBLE  PRECISION  CC < NM ) *  COL  * XNORM 
HAVE=H 
DO  2  1=1 fN 

DO  1  J=1 r N 
B  < I » J  >  =0 ♦ DO 
C ( I >  J ) =0 ♦ DO 
PP ( I f J ) =0 . DO 

AA(I»J)=DBLE<REAL<A(I»J>) > 

t  CONTINUE 

2  CONTINUE 
DO  17  1  =  1  *N 
COL-O* DO 
DO  16  J  =  1 » N 

C0L=DMAX1 <  COL » DABS ( DBLE <  REAL ( A ( I ?  J  )  ) )  )  ) 

16  CONTINUE 
CC  ( I ) =COL 

17  CONTINUE 
XNORM=CC  < 1 ) 

DO  18  I  =  1 >  N 

IF ( XNORM  . GT .  CC(I))  GO  TO  18 
XNORM=CC ( I ) 

18  CONTINUE 

*  COLUMN  NORM  IS  USED  TO  SEE  WHETHER  EXP(A)  NEEDS  REDUCTION  * 
M=0 

30  IF(XNORM*H  -  0.98)  3 t 20 >20 

*  EXP  <  A )  =  ( I- . 5#A )#*<-l)#(I+. 5#A )  * 

)|(  )|(  ]|[  ]||  ^  ^  ^  )|c  }|c 

3  DO  6  1  =  1  »N 
DO  5  J=1 >N 

DO  4  K  = 1 » N 

PP < I f  J ) =PP ( I r J ) +AA ( I » K ) *AA ( K * J > 

4  CONTINUE 
C  ( I »  J )  =~AA  ( I  r  J )  #H/2 . 0 

5  CONTINUE 
C< I » I )=C< I f I )+l .DO 

A  CONTINUE 

CALL  INVERT <  C  » N  >  B  ) 
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